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We present the first results for neutral kaon mixing using 2+1 flavors of domain wall fermions. A 
new approach is used to extrapolate to the physical up and down quark masses from our numerical 
studies with pion masses in the range 240 - 420 MeV; only SU(2)l X SU(2)r chiral symmetry is 

assumed and the kaon is not assumed to be light. Our main result is B% s (2 GeV) = 0.524(10)(28) 
where the first error is statistical and the second incorporates estimates for all systematic errors. 



The phenomena of CP violation is a central component 
of the Standard Model, in which CP violation is only pos- 
sible when all three of the quark doublets present in Na- 
ture interact. The Cabibbo-Kobayashi-Maskawa (CKM) 
flavor mixing matrix contains a single, physically mean- 
ingful phase which must describe all CP violating phe- 
nomena. For bottom mesons, it is possible to make a di- 
rect connection between the measured CP violation in B 
decays and this CKM phase. However, for K mesons, the 
system in which CP violation was originally observed, 
this connection is far more challenging. 

One begins with the measure of indirect CP viola- 
tion e K = 2.232 ± 0.007 ■ 10~ 3 Q, determined experi- 
mentally from the mixing between K° and K mesons. 
The operator product expansion relates ck to the QCD 
matrix element of a four quark operator Ovv+aa = 
{s^p,d){s^^d) + (s757 M ci)(s757 M d) between kaon states via 
a well known perturbative expression Q involving this 
CKM phase. This matrix element is parameterized by 
the renormalization scheme dependent parameter 

(K \Ovv +AA \k°) 
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Lattice QCD offers the only first-principles determination 
of Bk, which is essential to determine if the CP viola- 
tions observed in the B and K systems have a common, 
Standard Model origin. We describe a lattice QCD calcu- 
lation of Bk in which the most important errors present 
in earlier lattice results have been substantially reduced. 
We exploit the domain wall fermion (DWF) formulation 
with 2+1 dynamical flavors. This suppresses O(a) er- 



rors, both on- and off-shell, and also chiral symmetry 
breaking (measured by a (small) additive "residual mass" 
m rcs ). This allows us to renormalize Ovv±AA, multiplica- 
tively via a non-perturbative matching [113, H- Thus, 
we simultaneously avoid the complexity of lattice opera- 
tor mixing, avoid poorly convergent lattice perturbation 
theory and include the correct light flavor content. 

Alternative lattice approaches to Bk must treat either 
a chirality or taste mixing matrix and result in larger 
errors. For Wilson fermions a chirality mixing matrix 
can be determined using non-perturbative off-shell renor- 
malization, but large cancelations leave results impre- 
cise. While staggered fermion simulations successfully 
treat simpler quantities, current staggered results have a 
10-20% error for Bk due to large taste mixing [6j. 

By using large lattice volumes and meson masses as 
light as 243 MeV, we can determine the light quark limit 
with substantially improved accuracy. Instead of using 
a chiral perturbation theory (ChPT) which treats the K 
meson as light compared to the chiral scale, we evalu- 
ate the chiral limit using SU(2) L x SU(2) R ChPT and 
assume that only our pions are light. 

SIMULATION 

Our calculation is performed with a fixed lattice spac- 
ing and two space-time volumes, 16 3 x 32 and 24 3 x 64, 
using the Iwasaki gauge action 0] with j3 — 2.13 and 
the DWF action with a fifth dimension of size 16. Each 
ensemble uses the same dynamical strange quark mass 
am sea _ o.04 in lattice units. We use three 16 3 en- 
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sembles with degenerate up and down quarks of mass 
am]™ G {0.01,0.02,0.03} and two 24 3 ensembles with 
aTO sea g {0.005,0.01}. The ensembles, described in 



were generated using the RHMC algorithm [10| 
with trajectories of unit length. The 16 3 ensembles 
each contain 4000 trajectories, from which we omit 1000 
trajectories for thermalization. We perform measure- 
ments on 150 configurations separated by 20 trajecto- 
ries for each ensemble, calculating matrix elements of all 
possible pseudoscalar states with valence quark masses 
am val G {0.01,0.02,0.03,0.04,0.05}. We bin the data 
using up to 80 trajectories per bin to reduce the corre- 
lations between our samples. The am™ = 0.005 and 
0.01, 24 3 ensembles are composed of 4460 and 5020 tra- 
jectories respectively with measurements performed on 
the final 90 configurations separated by 40 trajectories. 
These data are binned into blocks of 2 configurations 
representing 80 trajectories. All possible pseudoscalar 
states are studied, composed of valence quark masses 



^ val 



e {0.001, 0.005, 0.01, 0.02, 0.03, 0.04}. 



We use the mass of the fl baryon, linearly extrapo- 
lated to mi = (m u + md)/2, and mx and m n treated in 
SU L (2)xSU R (2) ChPTto determine 1/a = 1.73(3) GeV, 
and unrenormalised masses am s + am rcs — 0.0375(16) 
and ami + am ICS — 0.00130(6). Since in our simulation, 
am s = 0.04 and am rcs = 0.00315(2), we must take into 
account our 15% too large input value of m s . We ob- 
tained the pseudoscalar decay constant af n — 0.0718(18) 
using Z^(P|j4o|0) = —iMpfp, and the axial current 
renormalization constant Za = 0.7162(2) (Note, the 
above errors are all statistical.) 

We use an established method [1, [I]], [l2} for evaluating 
the lattice matrix elements. Zero-momentum kaon states 
are created and annihilated using Coulomb-gauge-fixed 
wall sources at times 5 and 27 for the 16 3 volume and 5 
and 59 for 24 3 , with Ovv+aa inserted at each intervening 
point. Use of a combination of periodic and anti-periodic 
boundary conditions in time removes unwanted propaga- 
tion around the boundary, resulting in a good plateau for 
the ratio of matrix elements 



Slat 
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A sample fit is displayed in Fig. [TJ For each ensemble 
the pseudoscalar mass, Mp, decay constant, f P , and 
.B-parameter, Bp, are computed for each combination 
(m™ 1 ,™™ 1 ). We tabulate a portion of the 16 3 and 24 3 
results for B P and Mp in Tables HI and HT1 

Our lightest dynamical pion masses are 331 MeV. We 
must extrapolate our result for Bp to the physical value 
of mi — (ra u + md)/2, and we treat only the up and down 
quarks as light by using the SU(2)l x SU(2) r partially 
quenched ChPT (PQChPT) formula [H 
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FIG. 1: Typical plateau for the lattice .B-parameter for the 
pseudoscalar state made up from quarks of mass amj a = 
0.001, amf = 0.04, on the amf a = 0.005, amf a = 0.04, 24 3 
ensemble. 
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Bp 



= 0.01 

aMp 



Bp 



= 0.02 

aMp 



am? ea = 0.03 



Bp 



aMp 



1.1 
1.2 
2.2 

1,4 
2.4 
3.4 
4,4 



0.546(8) 
0.577(6) 
0.598(5) 
0.620(5) 
0.633(4) 
0.647(4) 
0.659(3) 



0.247(3) 
0.290(3) 
0.323(3) 
0.356(3) 
0.387(3) 
0.414(3) 
0.438(3) 



0.539(8) 
0.569(6) 
0.589(5) 
0.616(6) 
0.627(4) 
0.641(4) 
0.655(3) 



0.250(3) 
0.292(3) 
0.325(3) 
0.359(3) 
0.388(3) 
0.415(3) 
0.440(3) 



0.527(7) 
0.556(6) 
0.580(5) 
0.599(5) 
0.618(4) 
0.634(3) 
0.648(3) 



0.251(3) 
0.289(3) 
0.325(3) 
0.356(3) 
0.385(3) 
0.412(3) 
0.442(3) 



TABLE I: Bare pseudoscalar B-parameter Bp, and mass Mp, 
for the 16 3 volume. The fit range is 12-20 for Bp and 15-27 for 
Mp. In the first column, e.g., 1,2 denotes a meson composed 
of quarks with lattice masses (omJ al , am™ 1 ) = (0.01, 0.02). 

Here m x is the light valence quark mass, mi the light sea 
quark mass, A c h the chiral scale, Bp(Q), cq and c\ are 
m s -dependent low energy constants, k = Bomi / (Air f^) 2 , 
Bq is the constant in the expression m\ — 2Bom x , and 
we include m rcs in all masses entering these formulae, 
mjf is now the lowest scale that can dimensionally bal- 
ance higher order terms in the chiral expansion, and 
m\ jm\ will determine the suppression of successive or- 
ders. This is a better expansion parameter than the 
m 2 K /A 2 ch of SU(3) ChPT. We should note that the pseu- 
doscalar masses and decay constants are also well de- 
scribed by a SU(2) L x SU{2) R PQChPT analysis. 

Figure shows Bp versus m x together with the 
SU{2) L x SU{2) R partially quenched ChPT fit to the 24 3 
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0.001,0.001 


0.469(8) 


0.1402(9) 


0.470(5) 


0.1434(10) 


0.001,0.005 


0.491(7) 


0.1681(8) 


0.495(4) 


0.1707(9) 


0.005,0.005 


0.508(5) 


0.1916(8) 


0.512(3) 


0.1938(8) 


0.001,0.01 


0.514(6) 


0.1971(8) 


0.521(4) 


0.1995(9) 


0.005,0.01 


0.527(4) 


0.2172(8) 


0.531(3) 


0.2194(8) 


0.01, 0.01 


0.542(3) 


0.2400(7) 


0.546(3) 


0.2421(8) 


0.001,0.04 


0.601(6) 


0.3204(11) 


0.613(7) 


0.3234(11) 


0.005,0.04 


0.607(4) 


0.3329(8) 


0.611(3) 


0.3358(8) 


0.01, 0.04 


0.614(3) 


0.3482(7) 


0.616(2) 


0.3509(7) 



TABLE II: Bare pseudoscalar B-parameter Bp, and mass 
Mp, results for the 24 3 configurations. The fit range is 12-52 
for the B-parameter and 15-59 for the mass. 
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FIG. 2: Results for Bp together with the NLO partially 
quenched SU{2) L x SU{2) R ChPT fit to the 24 3 data plotted 
versus the light valence quark lattice mass am x . From top 
to bottom on the left-hand-side, the three curves are ami = 
0.01, 0.005 and am x respectively. The valence strange quark 
mass is fixed at its unitary value am v — am 3 = 0.04. While 
the statistical errors are large, the growing upward curvature 
in m x as the sea quark mass is increased from 0.005 to 0.01 
predicted by ChPT is visible. Some m x values are slightly 
shifted for clarity. 

data. The fit does not include correlations so the result- 
ing x 2 /dof = 0.14 is not a meaningful indication of good- 
ness of fit. The 16 3 unitary data, also shown in Fig. [2l 
are well described by a straight line which, if simply ex- 
trapolated to the physical limit gives a result about 6% 
larger than the more accurate chiral extrapolation that 
is possible if the smaller masses in the 24 3 simulation are 
used [13]. This SU(2) L x SU(2) R chiral extrapolation 
gives Bk = 0.565(10) for physical mj. (Bk is deter- 
mined at the correct valence value of am s = 0.0343 by 
linearly interpolating between am y = 0.03 and 0.04.) 
We have previously demonstrated 
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wrong-chirality mixing in our simulation is sufficiently 
suppressed that we may ignore it. Thus, we consider only 
multiplicative renormalization of Oyv+AA, and use the 
RI-MOM non-perturbative renormalization technique to 
match our lattice scheme to the MS scheme via B^ s = 



( Z Ovv + aa/ Z a) B K 



, 72 \ nlat 

SVV+AA I 

The technique performs well since domain wall 
fcrmions are off-shell improved. We evaluate the am- 
putated, four- leg and two-leg vertex functions Ae> vv+AA , 
A^4, and Ay in Landau gauge. Because of the relatively 
low lattice cut-off, A^ and Ay differ by 2%. We use their 
average and add their difference to the systematic error. 

We both quote lattice results in the RI-MOM scheme 
without perturbative error and convert to other schemes 
using the continuum NLO result [H, [l3]- We obtain 
Z%\{2 GeV) = 0.910(5)(13). The first error in paren- 
thesis is statistical and the second systematic. As shown 
in Figure [3l we use the RGI scheme as a scale invariant 
intermediate step to reveal and remove possible (ap) 2 er- 
rors. Only weak scale dependence is seen in the window 
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FIG. 3: A plot of Zg C ^, I (p 2 ) showing that the perturbative run- 
ning, removed from Zg^, 1 , accounts for most p 2 dependence. 



1.0 < (ap) 2 < 2.5, implying artefact-free perturbative 
behavior. We obtain Z^ 1 = 1.275(10)(25) by linearly 
extrapolating to (ap) 2 = to remove 0(a 2 ) effects. 

Conversion to MS is a 2% effect at NLO, consistent 
with 0(1) x a s /47T. While the error estimate could be as 
low as O(few) x (a s /4n) 2 , we add in quadrature the size 
of the NLO correction itself as a perturbative systematic, 
giving Z]^(2 GeV) = 0.928 (5) (23). For comparison, 1- 

loop lattice perturbation theory [l8| gives Z^ s (2 GeV) = 
1.007 with the difference likely due to slow convergence 
of lattice perturbation theory. 



RESULTS AND CONCLUSIONS 

We combine the bare Bk given above with these Z 
factors to obtain physically normalized results from 2+1 
flavor DWF at a" 1 = 1.73(3) GeV on a 24 3 x 64 volume, 
where the second error is the renormalisation system- 
atic: B K l {2 GeV) = 0.514(10)(7), and B^(2 GeV) = 
0.524(10)(13). The RI result involves no use of pertur- 
bation theory and has a reduced systematic error. 

We plot the MS result in relation to those quenched 1 1 
and two flavor 12 DWF results that could allow de- 
pendence on flavor content to be seen. (Of course dif- 
ferences arising from different lattice spacing and gauge 
action are also to be expected.) In addition, we include 
the 2+1 flavor staggered fermion result^ in Fig. [4] 

Finite volume chiral perturbation theory [2(| suggests 
that finite volume effects are negligible for all masses and 
volumes in our simulation except for am x = 0.001 where 
the effect may be 2%. However, since our fit is insensitive 
to this point, we adopt our 1%, 16 3 - 24 3 , ami = 0.01 
difference as an estimate of the finite volume error. 

Finite lattice spacing errors are likely larger, but also 
difficult for us to estimate. We can make use of the 
quenched, perturbatively renormalized results of CP- 
PACS[ig| also obtained for the Iwasaki and DWF lat- 
tice action which suggest a (poorly determined) scaling 
violation of size 3.5% at our coarser lattice spacing. We 
choose a 4% systematic error as the most likely estimate 
for 0(a 2 ) effects in our 2+1 flavor result. A second ap- 
proach to estimating discretization errors is to compare a 
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variety of presumably reliable quantities computed from 
our l/a — 1.73(3) GeV ensembles with experiment. For 
example, we find f„ and about 4% below experiment, 
a discrepancy consistent with our 4% error estimate. 

While we have interpolated to the physical valence 
strange quark mass, we have results for only one strange 
sea quark mass. We estimate the error resulting from 
our 15% too large m s by observing that for fixed valence 
masses (0.01,0.04) Bp increases by 3% when the light sea 
quark mass is changed from 0.03 to 0.02. Scaling this to 
the 0.0057 change needed for a single flavor of sea quark 
implies a 1% error. Finally we add a 2% chiral extrap- 
olation error by estimating the effects of NNLO ChPT 
as the 6% difference between the linear and NLO chiral 
limits in Fig. [2] scaled by mi/m s = 0.4 for am; = 0.01. 

Thus, we take our central value, which removes all 
quenching systematics, and add the 1% finite volume, 
4% scaling, 1% m s extrapolation, 2% ChPT and renor- 
malization error estimates in quadrature and obtain: 

B^(2 GeV) = 0.514(10)(25), (3) 
BjF(2 GeV) = 0.524(10)(28), (4) 
B K = 0.720(13)(37), (5) 

where the first error is statistical and the second is the 
estimated systematic. A recent review, including all lat- 
tice data then available, quoted a continuum limit value 
of Sf^(2 GeV) = 0.58(3) (6) [H]]. Our result is consistent 
with this and reduces both types of error substantially. 
This improvement arises because, using QCDOC com- 
puters, we have for the first time combined the correct 
dynamical flavor content with a lattice formulation with 
good chiral symmetry, O(a) improvement, control over 
operator mixing, non-perturbative renormalization and 
a new use of SU(2) L x SU(2) R ChPT for kaons. 

In order to reduce the significant discretisation system- 
atic error in our result for Bk , the RBC and UKQCD col- 
laborations are now doing simulations at a smaller lattice 
spacing, which will give quantitative control of this effect. 
A continuum NNLO perturbative calculation is required 
to convert lattice results to MS with better precision. 
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FIG. 4: We compare our 2+1 flavor results with earlier 
quenched[l{| and 2 flavor DWF 12] as well as 2+1 flavor stag- 
gered calculations @]. The quenched Iwasaki points show sta- 
tistical errors only while our point and the staggered point 
include renormalisation systematics. While our point lies be- 
low these Iwasaki results due to our improved chiral limit and 
flavor content, we expect similar a 2 dependence. 
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